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Abstract 



^ -Predicting the rate of escape and thermal structure of Pluto's upper atmosphere in preparation for the New Horizons 
Spacecraft encounter in 2015 is important for planning and interpreting the expected measurements. Having a moderate 
q Jeans parameter Pluto's atmosphere does not fit the classic definition of Jeans escape for light species escaping from 
^ "the terrestrial planets, nor does it fit the hydrodynamic outflow from comets and certain exoplanets. It has been pro- 
^ posed for some time that Pluto lies in the region of slow-hydrodynamic escape. Using a hybrid fluid/molecular-kinetic 
'model, we previously demonstrated the typical implementation of this model fails to correctly describe the appropriate 
temperature structure for the upper atmosphere for solar minimum conditions. Here we used a time-dependent solver to 
p | allow us to extend those simulations to higher heating rates and we examined fluid models in which Jeans-like escape 
pr"] expressions are used for the upper boundary conditions. We compare these to our hybrid simulations of the atmosphere 
y-H ^under heating conditions roughly representative of solar minimum and medium conditions as these bracket conditions 
Q_i expected during the New Horizon encounter. Although we find escape rates comparable to those previously estimated 
q .by the slow-hydrodynamic escape model, and roughly consistent with energy limited escape, our model produces a 
£j much more extended atmosphere with higher temperatures roughly consistent with recent observations of CO. Such an 
^ 'extended atmosphere will be affected by Charon and will affect Pluto's interaction with the solar wind at the New Hori- 
i zon encounter. Since we showed earlier that such models can be scaled, these results have implications for modeling 
exoplanet atmospheres for which the energy limited escape approximation is often used. 
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Recent spacecraft exploration of the planets and 
moons in our solar system and the rapid increase in the 
discovery of exoplanets has increased interest in atmo- 
spheric escape from planetary bodies. The Cassini space- 
craft is currently improving our understanding of atmo- 
spheric escape from Titan, while in 2015 the New Hori- 
zons (NH) spacecraft will flyby Pluto performing occulta- 
tion observations of its and Charon's atmosphere, and the 
Maven mission will orbit Mars studying the composition 
of its escaping atmosphere. Furthermore, newly discov- 
ered exoplanets with radii ranging from just a few times 
that of Earth to of the order of Jupiter have been discov- 



ered, and their atmosphere's also mod eled (La mmer et al. 



20091 : iMurrav-Clav et all 120091 : lYelleL I2004) . Here we 
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carry out simulations of Pluto's upper atmosphere that 
include atmospheric loss by thermal escape that can be 
tested against data to be obtained during the NH en- 
counter. By accurately describing the present loss rates, 
one can in principle learn about the evolution of Pluto's 
atmosphere. In addition, doing this accurately for a planet 
for which we will have in situ spacecraft data can, by scal- 
ing, guide our ability to model exoplanet atmospheres for 
which there will only be remote sensing data. 

The previous models of atmospheric escape for Pluto 
used the concept o f hydrodynam ic escape by adapting the 
critical solution in Parker ( 1964). Parker described the ex- 
panding stellar corona and stellar wind assuming the tem- 
perature and pressure go to zero at infinity, and showed 
that for this to happen the bulk velocity must increase past 
the isothermal speed of sound at a critical point that is de- 
pendent on the temperature and gravity. 

This model was subsequently adapted for pla netary at- 
mospheres to include heating by solar radiation (IMcNutt , 



ll989l:lHunten and Watsonlll982l) and was applied to Pluto 
dKrasnopolskvl . 1 19991 : Istrobell 12008b and Titan dstrobell 
2009). It is often referred to as the slow-hydrodynamic 
escape (SHE) model. The model required solving the 
fluid equations out to very large distances from the planet 
to enforce the necessary boundary conditions. However, 
it is known that at some finite distance the equations of 
fluid dynamics fail to describe the fl ow of mass, momen- 



tum and energy in an atmosphere (IVolkov et all 12011 



Johnson . 2010). The region of validity of the fluid equa 



tion is often classified using the Knudsen number, Kn = 
€/H, the ratio of the mean free path of the molecules, 
i = (V^trn) -1 , to the density scale height, H = G ^ /f2 , 
of the atmosphere. The fluid equations apply well where 
Kn <s. 1 so that many collisions occur over relevant length 
scales keeping the gas in thermal equilibrium. 

The concept of energy limited escape, which we 
will also examine here, is heavily used in m odeling es 



cape from early terrestrial planet atmospheres (ITian et al 



2009b a nd the growing body data on exoplanet atmo 



spheres ( Lammer et al. . 2009 : Valencia et al. . 201fJ) . For 
an atmosphere heated by EUV radiation, assuming that 
kinetic and thermal energy terms are small compared 
to gravity, and the thermal conduction is inefficient, the 
molecular loss rate is often approximated as 



Q 



U(r p ) 



(1) 



Here Q is the EUV energy supplied in the upper atmo- 
sphere, and U(r) = GMm/r is the gravitational energy. 
Without doing a detailed heating model one can use Q = 
nr EUV r l^ EUV ' wnere r EUV is the mean radius at which the 
radiation is absorbed and r\ and Few are tne heating ef- 
ficiency and solar energy flux respectively. This depends 
critically on tew which is sometimes assumed to be close 
to the visual radius so that r p can be used. The radius used 
to evaluate the gravitational potential is uncertain, with r p 
being a good option for terrestrial planets. 

The alternative to organized outflow is evaporative es- 
cape, by which the atmosphere loses gas on a molecule- 
by-molecule basis driven by conductive heat flow from 
below. The standard analytic model was originally devel- 
oped in Ijeansl d 1925b and is referred to as Jeans escape. 
The escape rates are found to depend predominately on 
the Jeans parameter A = GMm/rkbT, the ratio of grav- 
itational to thermal energy of the molecules. For large 
Jeans parameter the escape rate is obtained by integrat- 
ing the velocity distribution over the portion of molecules 
that are upward moving with speed in excess of the escape 



velocity. Assuming a Maxwell-Boltzmann velocity distri- 
bution the molecular escape rate cp and energy escape rate 
<pE from a altitude r are given by 



1 ? 
((f,) j = -nv • (1 + A) exp(-/l) • Anr 

<**>' = <'>''M£-§ + TTa 



(2) 
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Here v = s/SkhT/mn is the med molecular speed. To 
be consistent with the fluid equations, we also consider 
expressions which apply a drifting Maxwell-Boltzmann 
to include the bulk velocity u in the velocity dis tribution 
dYelleLl2004l : lTian et all Eoffll : IVolkov etailboilh . Typi- 
cally these equations are applied at a level called the exo- 
base where Kn « 1 , as there are few collisions above this 
level to inhibit a molecule from escaping. 

Accurate calculations of thermal escape can be ob- 
tained, in principle, by solving the Boltzmann's equation 
or, more often, by using th e Direct-Simulation Monte - 
Car lo (DSMC) model a s in iTucker and Johnson! d2009b 
and IVolkov et all (1201 lh . In this method, representative 
particles are used along with a model of molecular col- 
lisions to calculate the temperature, density and other 
macroscopic values from the velocity distribution in both 
the dense and rarified regions. In this respect, DSMC can 
model the exospheres with better accuracy than hydrody- 
namics model, but can only model the lower, denser re- 
gions of the atmospher e at great computati onal cost. Us- 
ing the DMSC method IVolkov et all d201 lb demonstrated 
earlier that for a monatomic or diatomic gas in the ab- 
sence of heating above some lower boundary for A > 3 
the escape is Jeans-like at the exobase (i.e. (/>/{<f>)j « 1.5 
for the range of Jeans parameters studied). For our model 
of Pluto, the lower boundary of the simulated domain has 
A = 22.8, and the exobase values exceed 4. Using the 
above results as guidance, we model the principal compo- 
nent in Pluto's atmosphere to obtain a description relevant 
to the NH encounter to test energy-limited escape, and 
to better understand the transition from Jeans to hydrody- 
namic escape. 

2. Model Description 

The steady state equations of mass, momentum, and 
energy have been used to study hydrodynamic escape. 
They can be solved to give the radial dependence of num- 
ber density n, outward bulk velocity u, and temperature T 
in the region in which th e atmosphere is col l isional. How- 



ever, as we have shown (Tucker et al., 2012; Volkov et al. 



20111) they cannot be used by themselves to determine the 



escape rate unless the Jeans parameter is very small in 
which case the solutions are somewhat insensitive to the 
boundary condition at infinity. The equations, neglecting 
viscosity can be written as 



d (\ 2 
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' dr r 2 
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= Anr q(r) 
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Here <p is the molecular escape rate through the one di- 
mensional atmosphere, p is the pressure (related via the 
equation of state, p = nki,T), m is the mass of N2 
molecule, G is Newton's gravitational constant, M is the 
mass of Pluto, is Boltzmann's constant, k(T) is the con- 
ductivity, C p - ^kb is the specific heat at constant pres- 
sure, and q(r) is the net heating/cooling rate per unit vol- 
ume. In Eq. ©, the first term on the left side is work done 
by adiabatic expansion and the second is the work done by 
conduction. 

For the conductivity, we use the power law k{T) = 
kqT s to approximate the temperature dependence. Some 
authors use an empirical fit for the conductivi ty, e.g., 
iHunten and Watsonl (1198211 used s = 1.12. While iMcNutt 
(1989) used 5 = 1 since it simplifies the analytic solu- 
tions to Eq. (©. In this paper we us e Kg = 9 .37 x 10~ 5 
J/m K and s = 1 to compare with ISfrobell ((2008) and 



because this is consistent with the variable-hard-sphere 
model for nitrogen molecule collisions, including trans- 
lation/rotatio nal coupling, used in DSMC models of the 
exosphere dTucker et all 120 1 2l : R^lkov et allboill) . 

The radiative heating and cooling model is adapted 
from lstrobell d2008h . using the same fixed mixing ratios of 
Pn 2 - 0.97, pcH 4 = 0.03, and pco = 0.00046 for comput- 
ing the heating and cooling rates. The energy fluxes and 
effective cross section for FUV ab sorption by CH4 and 
EUV absorption by N2 are given in Krasnopolsky (1999) 
for the various levels of solar activity at Pluto's perihelion 
(i.e. 30 AU). We obtain globally averaged heating rates 
for the FUV/EUV by applying Lambert-Beer's Law for 
a plane-parallel atmosphere separately for each species s: 
F s (r) = ^f*Ff exp(-r s (r)///), where F s (r) is the energy 
flux, r s (r) = f°° cr s n s (r)dr is the vertical optical depth, 
cr s is the cross section of species s, n s is the number den- 
sity of species s, and p = cos(60°) = 0.5 is used to ap- 
proximate spherically averaged heating (iStrobel and ZhuL 
1994 lApruzesa. 1 1980 r). Then the heating rate for each 



species is given by 

dF s cr s n s (r) 

q s (r) = —r~ - F s (r) 

dr U 



(7) 



Further heating and cooling is provided by methane near- 
IR absorption and CO rotational line emission both of 
which are non-LTE proce sses. Pa r amete rizations of these 
processes are given m IStrobell d2008h . fit to the de- 
tailed radiative-transfer model of Pluto's lower atmo- 



sphere |S_ttobel_andZhu ( 1996). The methane heating rate 

exp(-^y^) ergs per sec, 



per molecule is given as 8 x 10 

GMm 
rhTo 



where A - and = A(ro). The CO rotational cool- 
ing rate per molecule is given as 8 x 10~ 20 (-^) 2 ' 45 ergs 
per sec. 

These calculations assume a fixed efficiency for the 
various processes, yet we know that this assumption will 
break down when the gas becomes tenuous. To accom- 
modate this we assume that above Kn - 0.1 the heating 
is zero (i.e. efficiency drops to 0). Further, with the heat- 
ing per molecule constant, including heating up to the top 
boundary can cause numerical instability even though the 
heating rate per unit volume is small. In describing the 
heating, a cut -off is o f ten us ed where the tangential line of 
sight is unity Strobel ( 20081) . which in our case lies below 
the Kn = 0.1 level. 

Previously we directly solved the time-independent 
equations, iteratively solving them along with the heating 
and the DSMC e scape simulations until a consist ent so- 
lution is found dTucker et all 120121 : Istrobei 120081) . Here 
we reintroduce the time-dependence into the energy equa- 
tion and time-step until convergence. This is found to be 
more numerically stable for the boundary conditions that 
we w ill use here. We make the substitution £ = T 
as m dStrobel and ZhuL Il996l : IZalucha et all l2010h in the 
Eq. © to obtain a linear time-dependent PDE: 



d£ (s + l)T s ( cf> 
(it ntp \ 

k T s 2 cf> 
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nC p dr 2 



This equation is solved using the finite differences and 
the implicit, time-stepping scheme described in Appendix 
A. This second-order equation needs two boundary con- 
ditions, the first is the fixed lower boundary temperature 
T(ro) = SS.2K. The second is a restriction on the energy 
flux leaving the atmosphere 



W J3L d l 

s + 1 dr 



r top 
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(9) 
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This condition is derived by integrating Eq. © from the 
exobase to the upper boundary, and using <pE for the en- 
ergy flux at the exobase. 

The time stepping advances the temperature alone, 
next we find density, n, and bulk velocity, u, at the new 
time step. First u is updated using Eq. © and the previ- 
ous n. This new value of u is used to update n using the 
momentum equation © expressed as 



n(r) = n exp |- 



" £(hT + 



W _GMm ) \ 

— — dr\ 



k b T 



(10) 



As these two equations are coupled through n and u, we it- 
eratively solve them together to find a consistent solution. 
A few iterations are sufficient to get the fully consistent 
solution for T, n, and u. With these profiles determined, 
the heating components are recalculated. 

We proceed with the time stepping in two ways. As 



m lTucker et al.1 (120121) . we use the previous DSMC values 
of cp and <pE and time-step until the solutions are consis- 
tent with these values. Then we run a new DSMC sim- 
ulation from values obtained where the Knudsen number 
Kn = 0.1. This is far enough below the exobase, in the 
collision-dominated regime, to act as a lower boundary 
for the single component DSMC simulation, and above 
which heating can be ignored. The DSMC results for es- 
cape <p and <p£ are then used in the next fluid solution. 
After a few exchanges between the fluid and DSMC mod- 
els, a consistent solution is found for each specific heating 
case. Although non-thermal processes in the exobase re- 
gion contribute to and can dominate the atmospheric loss 
rate, we ignore these here, consistent with previous work 
on escape from Pluto. Such processes will be included in 
subsequent work. 

Because the DSMC calculations can be time consum- 
ing, for comparison, we also approximate the escape pro- 
cess as Jeans-like. That is, in each iteration we find the 
exobase, Kn — \, and compute the Jeans escape values 
(including bulk velocity) for <p and <pE for the next time- 
step. A consistent solution is obtained by stepping in time. 
Since we have shown earlier using the molecular kinetic 
model that the Jeans rate consistently underestimates the 
escape rate for the range of A of interest here, we also scale 
the Jeans escape rates by a value ~2 consistent with our 
earlier DSMC results. 

The Jeans escape boundary condition have been 
used in modeling planetary a tmo spheres, notably by 
Chassefierel d 1 996^ . I Yellel d2004f) and lTian et allfcpQSh . In 
modelin g the resp onse of Earth's thermosphere to EUV 
heating, iTian et all used the Jeans escape expression to 



determine the velocity at the upper boundary with the up- 
per boundary being the exobase, but fixed the heat gra- 
dient at the upper boundary for the energy equation. In 
mode ling extra-solar gas giants at small orbital distance, 
lYellel also used Jeans escape, with a drifting Maxwell- 
Boltzmann velocity distribution, to get an upper boundary 
velocity. But lYellel applied this at a fixed upper boundary 
of 3 planetary radii rather than at the exobase. 

To aid in convergence to steady-state, we estimate the 
heating q at the advanced time step using a Taylor expan- 
sion in £ as q{r,g +1 ) - q(r,t?) + J(f' +1 - ). This can 
help dissipate numerical oscillations between time steps, 
speeding up convergence. In calculating the derivative for 
the FUV and EUV heating we assume ^ = and 4^ = j 
as the optical depth t « crnH with H depending linearly 
on T. 

The lower boundary is set at 1450 km, consist ent with 
the oc cultation results and with the assumptions in Strobel 
( 2008h . Since we showed earlier that the atmosphere is 
highly extended, the upper boundary is set to 16000 km 
to accommodate the solar maximum case. For the radial 
grid, we need a fine resolution for the first 3000 km where 
the heating occurs, but a coarse grid is sufficient above 
this. Therefore, we use a grid that is equally spaced in 
1/r. With 1600 points in r we get a lower boundary dr of 
order 1 km and an upper boundary dr of order 100 km. 

3. Results 

3.1. Fluid-DSMC 

With the time-dependent solver, which is an improve 
ment on the steady-state iterations used in iTucker et al. 
(120121) . we are able to model higher solar heating rates. 
We present the two solutions found in that paper for 
no heating and solar minimum heating along with a 
new fluid-DSMC hybrid result for solar medium heat- 
ing in Figure [T] This allows us to bracket the expected 
heating rates at the NH encounter. The hybrid solu- 
tion includes rotational/ translational energy exchange 
consistent with th e assumed thermal conductivity (see 
Tucker et all d2012h ). 

Looking at Figure |2j a few qualitative differences be- 
tween the hydrodynamic solutions and our fluid-DMSC 
solutions are seen. Although the two types of calcula- 
tions give similar escape rates, they p roduce different tern 



perature pro f iles a lso discussed in (Vol kov et all 12011 



Tucker et all 120121) . In the DSMC the parallel and per- 
pendicular temperature differ above Kn « 0.1. More im- 
portant for the NH encounter, the exobase altitude for 
the fluid-DSMC is much higher suggesting a significantly 
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Figure 1 : Temp erature (right) and d ensity (left) profiles for the hybrid fluid-DSMC model of Pluto, for the special cases of no heating and solar 
minimum (from lTucker et al.l ( l2012h ). and our new result for solar medium. The temperate profiles are plotted with the scale labeled below; the 
density profiles are plotted on a logarithmic scale, with scale labeled above. The exobase altitudes are marked for each case. 
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Figure 2: Temper ature (right) an d density (left) profiles for solar medium-like condition, comparing the fluid-DSMC solution of this study with 
that of SHE from Strobell (2008), along with the fluid- Jeans solution describe in Section [X2l Marked on the fluid-DSMC solution is the exobase 
(r x ) as well as the transition between the fluid model and DSMC (r,, where Kn(r t ) = 0.1). Marked on the SHE model are the exobase (r x ) and the 
sonic point (r s ). 
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more extended atmosphere for Pluto. Further, the hy- 
drodynamic methods can produce an atmosphere that ex- 
pands faster than the sonic speed. However, we find the 
expansion is subsonic with the bulk flow being less than 
3% the speed of sound at the exobase. Further out, the 
bulk speed increases to approach a nearly constant value 
well below the speed of sound, consistent with a Jeans- 
like model. 

In the solar medium case, the transition from the fluid 
to DSMC domain is near 3000 km. In the region 0.1 < 
Kn < 1 .0 the agreemen t between the fluid a nd the DSMC 
is better than 5% as in iTucker et all d2012h . The DSMC 
solution begins to deviate from the fluid solution near the 
altitude at which the kinetic and internal energies of the 
molecules begin to separate, which is associated with the 
transition to non-equilibrium flow. It is important to note 
that this occurs below the exobase, so that the region of 
validity of the fluid equations should be terminated below 
the exobase. Our choice of Kn = 0.1 for connecting the 
DSMC solution to the fluid solution is at an altitude below 
which this separation occurs. 



Strobell ([2008) solutions result in exobase levels much 
lower than that found in our model. The correspond- 
ing exobase altitudes are ~0.56 - 0.7 of those calculated 
with the hybrid model. With these lower altitudes and the 
lower temperatures, the Jeans parameter at the exobase, 
A exo , are ~2.6 - 3.3 times larger than found in this study. 
Using these high A exo in Eq. © the Jeans escape rate is 
significantly lower than the escape rate found; hence the 
earlier conclusion that escape rates were orders of mag- 
nitude larger than the Jeans rate. In contrast the fluid- 
DSMC model results in an escape rate that is only mod- 
estly enhanced relative to the calculated Jeans rate, the ra- 
tio <pl(<p) j being 1.6, 2.0 and 2.3 for the cases considered. 
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Table 1: Fluid-DSMC results. E esc = if> E /<p is the average energy car- 
ried off per molecule; r exo is given as a ratio to r p = 1150 km. (cf>)j 
is computed using exobase values from the solution, <p L = Q/U(r) is 
computed using Q from the solution and r = 1450 km. 

The results in Figure [T] and Table Q] correspond to a 



net solar heating/cooling Q = £ x Anr 2 q{r)dr of 0.0, 3.8, 

and 7.8 xlO 14 ergs per sec. These can b e compared to 0.0, 
3.4, and 6.9 x 10 14 ergs per sec found in lStrobell d2008h . To 
give context for these heating rates, the upward flow of en- 
ergy across Pluto's 1450 km level is greater than 10 20 ergs 
per sec if it were escaping hydrodynamically (i.e. A — > 0). 

Table Q] gives many of the relevant values o f the fluid- 
DSMC hybrid solution in iTucker et aD d2012h including 
the new solar medium-like case calculated here using the 
time dependent model iteratively with a DSMC calcula- 
tions. It is seen that the escape fluxes are comparable to 
the Jeans flux and are also, quite remarkably, comparable 
to our estimate of the energy limited rate with the excep- 
tion, of course, of the Q = case. This result alone is 
important, as it is often presumed that the Jeans rate and 
the energy-limited rate are two limiting cases for atmo- 
spheric lOSS. 



As in ITucker et al.l (I2012h . for no heating our escape 
rate is a order of magnitude lower than the SHE solu- 
tion. For solar minimum we are 20% below the SHE so- 
lution, and for solar medium we are slightly above. In 
ea ch of th e se cas es we deposit somewhat more heat than 
in IStrobell (I2008h due to differences in atmosphere struc- 
ture. Scaled to the same heating rates the fluid-DSMC 
escape rate would be somewhat smaller than the hydrody- 
namic model as the DSMC produces higher temperatures 
in the exosphere region. 

3.2. Fluid-Jeans 

The fluid-DSMC solutions demonstrate that the es- 
cape rate and energy escape rate are within a factor of the 
order of unity of the Jeans rate calculated from Eqs. (0 
and (O using the temperature at Kn = 1. Therefore, in 
Figure |2] and Table |2] we present results using an upper 
boundary condition that is equal the Jeans escape and en- 
ergy loss rates, and then doubling those. It is seen that 
the fluid- Jeans solution is in many respects close to that 
of fluid-DSMC. Therefore, for Pluto's atmosphere under 
solar heating conditions likely occurring at the NH en- 
counter, approximating the escape process as Jeans-like 
reproduces many of the properties of the full solution that 
we are interested in. Although the fluid equations may not 
be exact up to the exobase, using them along with Jeans 
escape boundary conditions allows us to carry out param- 
eter studies that are impractical with the full fluid-DSMC 
model. 

From Table|2]we note that the escape rate is consistent 
with the fluid-DSMC model, but the other properties of 
the atmosphere differ, and this difference increases as the 
heating rate increases. The atmosphere is more extended 
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Figure 3: Temperature profile found in the fluid- Jeans model for the 4 special cases discussed: no heating, solar minimum, solar medium, and 
solar maximum. The top of each curve represents the exobase altitude. 



than in the fluid-DSMC model, allowing it to absorb more 
EUV heating. The Jeans parameter at the exobase is lower 
than that found in the fluid-DSMC model, which also can 
be attributed to the more extended atmosphere. 

In an attempt to more accurately model escape based 
on the DSMC we also used <p = 2 x (<p)j and (pE = 
2 x (<Pe)j for the boundary condition (note that the en- 
ergy per molecule, E esc , is unchanged by this scaling as 
found in the fluid-DSMC results). The most notable effect 
of this is to lower the exobase altitude and increase A exo . 
These results now more closely resemble the fluid-DSMC 
results in Table [TJ Surprisingly, ignoring the no heating 
case, the converged escape rate is barely affected by the 
enhancement of the Jeans boundary conditions. The same 
statement applies if we increase the enhancement to 3 or 4 
times Jeans. Since the escape rate is roughly given by the 
energy-limited approximation, this is consistent with the 
loss rates being insensitive to the size of the enhancement 
in the Jeans boundary conditions. Therefore, the other 
properties of the atmosphere adjust to ensure this. How- 
ever, we would not have had an estimate of the enhance- 
ment of the escape rate relative to Jeans in the absence 
of the fluid-DSMC simulations. Our early DSMC simu- 
lations indicated that for the case of no heating above the 
lower bo undary of the simu lation region, the enhancement 
is ~1.6 (iVoikov et all 1201 lh . How it increases with heat- 



ing remains an open problem. 

As we increase the heating and the atmosphere gets 
further extended, the Jeans parameter decreases to about 
3 at the exobase for the solar maximum case using the 
Jeans boundary co nditions. This is within the region that 



Volkov et al.l (201 1) found that marks the transition to hy- 



drodynamic escape for N2 at the relevant temperatures. In 
that study, it was shown that for the lower boundary at 
small Kn, the exobase moves out to infinity as the Jeans 
parameter approaches 3. Using twice the Jeans expres- 
sions in Eqs. © and (f3]> as the boundary conditions, the 
Jeans parameter in the upper atmosphere increases keep- 
ing the solution from entering the transonic regime. 

Based on the enhancement versus heating rate found 
in Table [T] we estimate the actual enhancement for solar 
maximum condition to be close to 3.0 times the Jeans rate. 
If we accept the assumption that the enhancement will be 
3.0 times Jeans for the solar maximum case then we can 
use the fluid- Jeans model to discuss the atmosphere. From 
Table[3]we observe the heating rate has about doubled rel- 
ative to solar medium, and the escape rate has slightly 
more than doubled. The average energy carried off per 
particle has decreased by ~30% relative to solar medium, 
but this is a property of the fluid-Jeans model that isn't 
representative of the fluid-DSMC model. The energy lim- 
ited escape rate 4>l - 5.61 x 10 27 for this Q, using U(tq), 



1 



Using lx Jeans 





none 


min 


med 


Q (10 14 ergs per sec) 





3.86 


7.87 


<p (10 27 per sec) 


0.035 


1.14 


2.58 


E esc (10~ 3 eV) 


15.5 


16.3 


14.6 


t"exo/fp 


3.41 


5.99 


8.10 


^■exo 


8.71 


4.72 


3.84 


Using 2x Jeans 








none 


min 


med 


Q (10 14 ergs per sec) 





3.83 


7.79 


<p (10 27 per sec) 


0.054 


1.17 


2.58 


E esc (10~ 3 eV) 


15.2 


15.4 


13.7 


^exo/^p 


3.34 


5.30 


6.84 


^■exo 


9.02 


5.55 


4.72 



Table 2: Results of fluid model using Jeans-like escape from the exo- 
base. First results using Jeans escape values for tp and <p E , while 
the second set uses twice the Jeans values to better capture the en- 
hancement of escape relative to Jeans. r em is given a ratio relative to 



1150 km. 



solar maximum 

lx 2x 



3x Jeans 



Q (10 14 ergs per sec) 


15.8 


15.7 


15.6 


<p (10 27 per sec) 


5.81 


5.80 


5.79 


E esc (10- 3 eV) 


11.7 


10.9 


10.4 


Fexo/fp 


12.4 


9.8 


8.9 


^■exo 


3.01 


3.95 


4.48 



Table 3: Results for fluid- Jeans model for solar maximum case. Given 
are results for lx, 2x, and 3x the Jeans value as a boundary condition. 
This last value is an estimate of the enhancement that would result 
from DSMC models, inferred from the results in Table 1. 



is close to our estimated escape rate. 

In Figure we plot the results for the four cases con- 
sidered here using the Jeans criteria in Eqs. (O and ([3]> at 
the upper boundary. The no heating case was initiated us- 
ing an isothermal profile and time stepped to steady state. 
Then we incrementally increased the solar flux rates to 
obtain the other 3 profiles in turn. As the exobase alti- 
tude increased with heating rate, the upper boundary al- 
titude had to be increased resulting in our use of 16000 
km, about 14r p , to model the solar maximum condition as 
discussed earlier. We see that unlike the exobase level, the 
altitude of the peak in the temperature does not increase 
significantly with heating rate, changing from 2124 km to 
2250 km from solar minimum to solar max with tempera- 
tures of 11 OK and 132K respectively. The heating peak 
is below the temperature maximum, and also does not 
change significantly with heating rate, going from 1767 



km to 1780 km. This heating peak is near where the opti- 
cal depth of the CH4 FUV is unity. The solutions become 
nearly isothermal at the exobase consis tent wi t h the fl uid- 
DSMC model, but different from thelParker! d 1964b and 
slow-hydrodynamic models JMcNuttLll989l : IWatson et all 



1 98 11 : iKrasnopolskvl . 1 1 9991 : IStrobelL 120081) that require the 
temperatures to decrease to zero. 

In his fluid dynamic studies of an early Earth atmo- 
sphere using the Jeans ra te for the upper boundary con- 
ditions, iTian et al.l (120081') found two regimes of escape 
in in response to increased UV heating. The transition 
from hydrostatic equilibrium to what he calls hydrody- 
namic escape coincides with adiabatic cooling becoming 
non-negligible when compared to conduction or radiative 
heating. In our simulations for Pluto, adiabatic cooling is 
always an important cooling mechanism in the upper at- 
mosphere. Further, we find similar changes in exobase al- 
titude, temperature and bulk velocity with increased heat- 
ing as in that paper's hydrodynamic escape regime, even 
though our escape is described as a kinetic process. That 
adiabatic cooling is always important in our simulations 
is likely due to the much smaller gravitational binding en- 
ergy and the similarities mentioned are likely due to his 
use of the Jeans escape rate at the upper boundary, which 
we showed here can be a reasonable approximation to the 
kinetic boundary conditions. 

In Figure |4] we compare the escape rate versus the 
net heating for the fluid-Jeans model to the other calcu- 
lations. The fluid-DSMC results are nearly the same, con- 
firming that the fluid- Jeans model can be used to estimate 
the escape rate for this range of Jeans parameter. With the 
exception of the no heating case, the SHE model results 
show a similar trend but are 43%, 16% and 7% higher than 
the fluid-Jeans for solar minimum, medium, and maxi- 
mum respectively, with the latter approaching blow-off. 

Also shown in Figure 0] is the energy-limited escape 
rate compared to our escape rate. The good agreement 
between the two suggests that most of the energy being 
deposited is going into removing the molecules from the 
gravitational well. The sources of cooling besides adia- 
batic expansion and loss to escape are conduction through 
the lower boundary and CO cooling (already included in 
Q), both of w hich are bounde d, as found in the hydrody- 
namic model (IStrobell. 120081) . The slope of the energy- 
limited escape rate can be better matched to our simula- 
tions by including the internal energy (i.e. C p (Tq - T exo )) 
in the denominator, and the vertical offset can be partly 
compensated for by subtracting the loss to conduction 
from the total energy deposited. 

We also compared estimates of the hydrodynamic 
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Figure 4: (Left) Escape rate t/> (solid line) plotted versus net h eat deposited Q from the fluid- Jeans model. Overlaid are escape rates from the fluid- 
DSMC model as black diamonds, and from the SHE model dStrobelL 2008) as x's. Also plotted is the energy limited escape rate <f> L = Q/U(r ) 
(dashed line). (Right) Comparing energy loss mechanisms versus net heat deposited from the fluid- Jeans model. 



escape rates in iHunten and Watson ( 1982) and iMcNuttl 
( 1989), and the energy-limited escape rate lLammer et all 
(2009) for the same solar flux. We calculated the heating 
by CH4 assuming a fixed mixing ratio of 3%, a globally 
averaged incoming solar energy flux, S = tjF, of 1x10" 3 
and 4x10" 3 ergs per cm 2 per sec (representing near so- 
lar minimum and maximum conditions respectively), an 
absorption cross-section for CH4 of 1.8xl0~ 17 cm 2 , and 
primarily nitrogen atmosphere with the values for conduc- 
tivity used in this paper. 

From Table |4] we see that neither of the simple 
hydrodynamic models produces a similar estimate to 
the fluid-Jeans model fo r either heating case. The 
Hunten and Watsonl (119821) model is an upper b ound, so 
it con sistently overestimates the escape rate. The lMcNutt 
(1989) is a good estimate for the smaller heating case, giv- 
ing a similar estimate of both escape rate and where the 
heating layer occur, but is quite different for the strong 
heating case. While these models can give good estimates 
of the escape rate for some cases, the solutions for the 
structure of the atmosphere are quite different from that 
found using our hybrid fluid-kinetic model. 

The energy-limited escape model is quite good if one 
knows the net heat deposited in the atmosphere, which de- 
pends on the atmospheric size. But without detailed mod- 
eling, one can only estimate the amount of energy lost 
through conduction and escape or deposited via heating. 
In Table |4] we consider using the energy-limited escape 
model with 4 estimates of the net heat deposited. First we 
use the net heating found from the fluid- Jeans model, then 
used Q = nrjS with the location of the heating layer be- 



ing the location of unity, parallel or vertical, optical depth 
(these altitude could also be estimated from observation), 
and finally with the heat deposit at the lower boundary as 
if no estimate were available. Knowing the net heating or 
the location of the heating layer provided a good estimate 
for the highly heating case, but for the near solar minimum 
case both overestimate the escape as a significant amount 
of heat is conducted through the lower boundary. Using 
the parallel optical depth provides a better estimate of the 
heat deposited. 

We investigated the sensitivity of the model to 
changes in some model parameters. If the atmosphere 
were a monatomic gas (C p = 5/2^), then E esc is about 
halved, primarily due to no internal energy being carried 
off, but the other properties are only significantly affected 
when the heating exceeds solar medium conditions. For 
solar medium, the escape rate increases by ~10% and the 
exobase is lower by a few hundred kilometers. In hydro- 
dynamic models it is often assumed that the energy carried 
off by escape is negligible and a positive non-zero value 
can only decrease the escape rate. If we use 4>e - (keep- 
ing <p as the Jeans value), then the escape rate increases by 
~5-10%. When the solar heating is above solar medium, 
corresponding to when (pE normally becomes a compara- 
ble energy loss process to conduction or CO cooling, the 
exobase altitude is lowered and the exobase temperature 
raised by as much as 20%. Therefore, even though ne- 
glecting (pE still produces a good estimate of escape rate, 
it should be included when interested in the structure of 
the atmosphere. 

Another consideration is our definition of Knudsen 
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s 


= 1 x Kr J erj 


*s cm 2 


sec 1 




(f>(l0 27 s" 1 ) 


rh/ro 


£(10 14 ergs s" 1 ) 


Fluid-Jeans 


1.06 


1.22 


4.68 


Watson 


15.4 


3.66 


35.4 


McNutt 


0.84 


1.17 


3.69 


riLc, (a) 


I. by 




4.00 


bLb (b) 


1.66 


1.32 


A ZTO 

4.62 


bLb (c) 


1 A 1 

1.41 


1.22 


i r\A 

3.94 


bLb (d) 


0.95 


1.0 


2.64 


S 


= 4 x 10" 3 erj 


*s cm" 2 


sec" 1 




0(1O 27 s" 1 ) 


rh/ro 


2(10 14 ergs s" 1 ) 


Fluid-Jeans 


8.46 


1.25 


23.3 


Watson 


32.1 


2.73 


76.8 


McNutt 


3.22 


1.28 


17.4 


bLb (a) 


8.40 




23.3 


bLb (b) 


7.30 


1.39 


20.3 


bLb (c) 


5.91 


1.25 


16.4 


bLb (d) 


3.81 


1.0 


10.6 



Table 4: Comp aring the fluid- Jeans model. lHunten and Watsonl l ll982h . 

and energy-limited escape (ELE). The values for each 
given are the escape rate t/>, the location of the heating layer in terms 
of the lower boundary r/,/ro (in the case of the fluid- Jeans model we 
give the location of the peak heating, which is also the unity verticle 
optical depth), and the net heating Q. The ELE model was applied in 
4 ways: (a) using Q from the fluid- Jeans model, (b) n, is altitude of 
unity parallel optical depth, (c) n, is r E uv from the fluid- Jeans model, 
(d) and r h is r as if no estimate were given. 



number used to determine the exobase. The relevant scale 
length in atmosphere is determined by the density gradi- 
ent as H n = -n/(dn/dr). In the above simulations we 
have used the isothermal approximation H = kt,T/mg(r). 
Using Kn - £/H n to determine the exobase results in a 
small correction (as H - H n though out the atmosphere) 
and has a negligible effect on the solution and escape rate. 
Since our atmosphere is highly extended one might con- 
sider Kn - £/r, which is appropriate for comets, where the 
density drops off rapidly as 1 /r 2 . Using this definition, the 
exobase altitude increases to 4.7, 9.3, 13.3, and 21.6 times 
r p for special cases. These values are about 50% larger as 
compared to those in Table [2] and [3] For solar maximum 
such an estimate would place the exobase beyond the orbit 
of Charon. 

3.3. Similar Escape Rates for Small Jeans Parameters 

This study has shown that the escape rates of the fluid- 
Jeans and fluid-DSMC are close to those found in the hy- 
drodynamic models when there is heating, and become 
closer as the heating increases. We posit that this can be 
explained by the Jeans parameter being lowered at the up- 



per boundary of the heating layer, and that for small Jeans 
parameter only a small range of escape rates lead to valid 
solutions as discussed below. 

In Figure [5] we plot the escape rate versus Jeans pa- 
rameter at the lower boundary for various models of es- 
cape. For these calculations, we rewrite the fluid equa- 
tions in a non-dimensional form, where the lower bound- 
ary conditions are given by the Knudsen number, Kno, the 
Jeans parameter, Aq, and the Mach number (scaled to the 



isother mal speed of sound), Mao, as s hown inlVolkov et al 



(120111) . In Figure [5] we compare the IParkerl (1 1964) solu- 
tion, which is the critical transonic solution that the many 
hydrodynamic models are based upon, to the fluid- Jeans 
so lution of this pap er. Also shown is the escape rate given 
by iGruzinovl (1201 11) in his analytic solution to these equa- 
tions. These were shown to match well with the tran- 
sition in escape rate found in the DSMC simulations of 



VolkovetalJ(l201lh . 

We see that the Parker and fluid- Jeans solutions be- 
gin to converge as the Jeans parameter decreases. As 
we increase the heating, the atmosphere expands and the 
Jeans parameter above the heating layer is decreasing, so 
the two models should produce similar escape rates. The 
fluid-DSMC solution lies between the Parker and fluid- 
Jeans as it is enhanced relative to Jeans but not a sonic 
solution. 

Solutions to the full hydrodynamic equations of 



Parkerl (119641) can be very sensitive to the escape rate, 
especially as the Jeans parameter decreases, and not all 
choices of escape rate can lead to valid solutions. The 
shaded region in Figure [5] gives the range of escape rates 
that support an exobase. The most important feature of 
this region is the change from containing static solutions, 
Mao - 0, to becoming narrow for small values of Jeans 
parameter. 

Both the fluid-Jeans and Parker's solution lie within 
this region, as does the slow-hydrodynamic escape solu- 
tion. Then for both solutions to have an exobase, they 
must have similar escape rates for small Jeans parame- 
ters. As stated above, even though the escape rates are 
within a factor of order unity, the upper atmosphere can 
have very different structures, from a super-sonic, rapidly 
cooling solution, to the near isothermal solution near the 
lower edge of the narrow region. 

When the atmosphere is heated this argument must 
be applied above the heating layer. Applying this at 
Kno - 0.1, where we have cut-off the heating, the nar- 
row region occurs for Ao < 6.7. We find that this value is 
achieved near solar medium heating. Therefore for solar 
medium and maximum heating the slow-hydrodynamic 
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Gruzinov 

Fluid-Jeans 

Parker 




Figure 5: Solutions from I Gruzinov! J201 ll> and lParkerl d 1964l> . compare to our fluid- Jeans solutions without heating for Kjiq = 10 5 (about Pluto's 
1450 km value). The shaded region is the range of escape rates that support an exobase. 



escape model and our hybrid models produce similar es- 
cape rates, even though the solutions are very different in 
describing the extended atmosphere. 

4. Conclusion 

We compared two methods for modeling the up- 
per thermosphere and characterizing esca pe. Having 



shown earlier that such models can be scaled Volkovetal 



(1201 lh . we considered Pluto's atmosphere as an example 
because of the anticipated NH encounter. We found that, 
in the presence of significant heating, the escape rate did 
not change significantly from that obtained in certain fluid 
models. However, many other characteristics of the atmo- 
sphere are affected, namely the exobase altitude and the 
existence of a sonic point. 

Since the fluid-DSMC model described earlier re- 
quires no assumption on upper boundary conditions, we 
used this as a base line for testing fluid models that re- 
quire assumptions about the particle and energy flux at 
the upper boundary. By iteration, the escape rate and en- 
ergy escape rate were shown to be continuous throughout 
the simulation domain as were the density, temperature 
and flow speed. We obtain a solution out to 10's of plane- 
tary radii. The escape rate found is approximately 2 times 
Jeans, but this value changes slowly with the amount of 



heat deposited in the upper atmosphere, as well as with the 
temperature at the lower boundary of the simulations, as 
shown earlier. The arrival of New Horizons should occur 
when the FUV/EUV flux is between the solar minimum 
and medium cases considered in this paper. 

The fluid- Jeans model uses the upper boundary con- 
ditions for Jeans escape as an approximation to use of the 
DSMC model in the exobase region. It was shown that 
the resulting escape rate is close to that given by the fluid- 
DSMC model. But we also show here that the small dif- 
ferences in escape rates do not imply that the temperature 
and density structures are similar. By scaling the Jeans es- 
cape rate by a factor between 2 and 3 the results become 
closer to the Fluid/DSMC. 

We also showed that the energy limited escape rate, 
often used for estimating atmospheric loss from exo- 
planets, was surprisingly accurate over a broad range 
of heating conditions. Since the details of the atmo- 
spheric structures of exoplanet atmospheres are usually 
not known, we also show that in applying this estimate one 
should use a radius associated with the FUV/EUV absorp- 
tion depth or, better, use that depth at which the parallel 
optical depth is unity. This depth can be roughly estimated 
by scaling to the visible radius using the ratio of the cross 
sections, the visible light scattering cross section to the 
FUV/EUV absorption cross, along with a scale height es- 
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timated using the equilibrium temperature. In this way the 
results here can be applied when modeling the evolution 
of exoplanet atmospheres. 

For the heating rates considered here we have not seen 
the transition from Jeans-like to hydro-like escape. How- 
ever, the resulting highly extended atmosphere for NH en- 
counter conditions means that the direction of the model- 
ing effort for Pluto needs to change. That is, in thermal 
escape calculations one terminates the energy deposition 
at very high altitudes as the heating efficiency of the back- 
ground gas due to solar radiation absorption drops rapidly 
with increasing altitude in the exobase region. However, 
the molecules in this the upper regions of this highly ex- 
tended atmosphere do absorb radiation. Although this 
might not efficiently heat' the atmosphere, a fraction of it 
does pr oduce non-thermal e scape, as described for other 



bodies (IJohnson et all 120081) . In addition, the interaction 



with the solar wind will now occur at much larger dis- 
tances from Pluto than had been pre dicted and can more 
efficie ntly strip the atmosphere (see iBagenal and Mcnutt 
For the extended atmosphere described here, 
these non-thermal processes at the NH encounter condi- 
tions could be comparable to thermal escape estimates and 
determine affect the NH observations. Such modeling is 
now in progress. 



would be in the case of equal spacing. The correct er- 
ror term is proportional to (Ar„, - Ar m _i )£"(?). While the 

quantity (Ar,„ - Ar,„_i) might be small compared to (Ar m ), 
the effect of this error in the model is to over-estimate the 
conductivity. To correct for this we use the finite differ- 
ences: 



$ V) - 2<- » - TZ bm+1 + \,-... \r " TT<Sm-\) 

LAt iti LA I n , 



2(Ar m -Ar m _i) , 
Ar,„Ar m _i 



Ar,„ 



Ar m _i Ar„, 



(A.l) 



f ( r ) — (t f=~^m+l a, A r T^ m + 

Ar,„Ar,-- tAr m ar m -i 



1 



'mAr m _i Ar„,_iAr„, 



(A.2) 



The error term on is now Ar„,Ar m _i£'"(r), now second 
order accurate and not affecting the conductivity. The er- 
ror term for f" is (Ar m - Ar m _i)£"'(r). While the finite 
difference approximation for f " is not second order accu- 
rate, this is the minimal error term one can get from the 3 
points used. It is more important that we capture conduc- 
tion correctly as to ensure a constant energy flux through 
the upper domain where there is no heating and the grid 
spacing is large. 

We adap t the implicit scheme proposed by 
Zalu cha et all 42010), using the 



d2010h . using 
solve the diffusion equation 



same notation, to 
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Appendix A. An implicit time stepping scheme with 
varying step size 

In this appendix, we describe the numerical equa- 
tions and techniques used to solve the diffusion equa- 
tion on an unequally spaced grid. We define the grid r m , 
m = l,...,M and a spacing Ar m = r m+ \ - r m . Define 
Ar„, - (r m+ i - r m -i)/2 as the average spacing at grid point 
r m . Last define % m as the finite difference approximation 
to £(r m ). 

First we describe approximate the derivatives d^jdr 
and d 2 £/dr 2 using unequal spacing. Using the points 
£ m +i, % m , % m -\ to approximate d^/dr, one might use the 
center difference scheme £'(r) « (£ m+ i -£ m _i)/(2Ar). But 
this scheme is not second order accurate, 0(Ar ), as it 



<9f d£ d 2 £ 



(A.3) 



on a domain r € (ro, r top ) along with the boundary condi- 
tions 



dr 



= F 



top 



top 



(A.4) 
(A.5) 



The coefficient B, C, G and source term 5 are functions of 
r and t. We make two modifications to the scheme. First 
we use our finite differences for unequal grid spacing, and 
second we modify the upper boundary condition to more 
accurately restrict the energy flow. 

Implementing the time grid ?, = iAt for i = 0, 1 , . . . 
and the notation %(r m , ?,) = g' m . Our implicit scheme is 
now obtained by substituting the finite differences, Eqs. 
(IA.1I ) and (IA.2I) . evaluated at the advanced time step into 
the diffusion equation. For the time derivative we use 
(£m _ £m)/A?- Solving for the advanced time step then 
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requires solving a tridiagonal system of equations 



(A.6) 



The lower boundary condition %\ - £o is then expressed 
as a\ = 0, Ai = 1, a>\ = 0, q\ - £o- The interior point (m 
= 2,. . . ,M-1) are expressed as 



Gtm — 



Ar m _iAr m I 2 



/im — 1 — A? I B,„ + 



1 



Ar m Ar m -\ 
((Ar m - Ar m _i)C m - 2G m ) 



-At 



Av m Av m 



Ar m -i 



q m = g m + AtS, 



(A.7) 



For the upper boundary we modify the procedure of in- 
terior points to enforce the condition on $ at the up- 
per boundary. We substitute Eq. (1A.5I ) in directly into 
Eq. (IA.3I) for the first derivative. Using the Eq. (IA.2I ) 
with m = M for the second derivative at the upper bound- 
ary introduces the phantom point £m+i> but by equating 
Eq. (IA.1I ) to Eq. (1A.5I ) we can remove it from the equa- 
tion. The result is expressed as 



a M = -AtG M 



Am - 1 - At ( Bm + Gm- 

oj m = 



1 + Ar M IAr M -\ 
Ar M Ar M -i 

2 - Ar M /Ar M -i 



Ar M Ar M 



<7m = i? M + At\S M + F,op\C M + 
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